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Recent theoretical advances in quantum electrodynam- 
ics (QED) proposed that arrays of coupled cavities, con- 
taining some form of optical nonlinearity, allow for a re- 
alization of a phase with strongly correlated excitations 
and photons [lj, |2|. The optical nonlinearity might be 
achieved for instance by coupling the cavity photons to a 
two-level system (TLS) in the form of an atom or an 
atom-like structure. Individual accessibility and high 
tunability of the parameters make coupled cavity sys- 
tems ideal candidates for quantum simulators and, prob- 
ably even more intriguingly, for applications in the field of 
quantum information processing. Initial theoretical work 
[3|— 15|] indicated a quantum phase transition of light from 
a localized Mott to a delocalized superfluid phase. Sub- 
sequently, the quantum phase transition [6|-|20| and the 
polaritonic nature of the constituent particles [ill, l20j 
have been investigated in great detail. In single-cavity 
QED experiments, an astounding high level of control 
has already been achieved, allowing the observation of 
numerous fascinating phenomena based on the nonlin- 
caritics at the single photon level [2l| - l37j . Importantly, 
in realistic experiments, these cavity QED systems are 
out of equilibrium, as they are driven by external lasers 
and are inherently susceptible to photon loss. Therefore 
the next step, having in view a scalable system of cou- 
pled cavities, is to understand the rich physics promised 
by a laser-driven and dissipative system of such coupled 
cavities, when out of equilibrium. 

Only very recently the theoretical investigation of such 
laser-driven and dissipative systems of coupled cavities 
has begun. The main research focus has been on the 
statistics of photons emitted from cavities with Kerr non- 
linearity, which are coherently driven by an external laser 
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FIG. 1. (Color online) Sketch of two laser-driven and dis- 
sipative coupled cavities (a). Each cavity confines photons 
and contains a two-level system. The coupling of the cavities 
arises due to the finite spatial overlap of their photonic wave 
functions as indicated by the yellow area. The first cavity is 
driven either incoherently (b) or coherently (c), with rate T + , 
see main text for details. The photon dissipation out of the 
second cavity occurs at rate F~ . 



source [38l - l45l ] . Steady-state entanglement of three co- 
herently driven and dissipative cavities has been studied 
as well [46[ . The time evolution of the population imbal- 
ance between two coherently driven cavities with Jaynes- 
Cummings (JC) nonlinearity |47| has been investigated 
in Ref. l48l In addition, the non-equilibrium quantum 
phase transition of a dissipative and coherently driven 
cavity system with Kerr nonlinearity has been analyzed 
using a mean field approach [49j. Finally, investigations, 
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within a slightly different context, on the non-e quili brium 
Bose-Hubbard model can be found in Refs. [5(M53l 

In this paper, we investigate excitations of laser-driven 
and dissipative coupled cavities with JC nonlincarity. To 
focus on the physical effects of interest, we introduce a 
'minimal model' consisting of two coupled cavities with 
a drive term applied to one two-level system and loss via 
photon emission from the other, see sketch in Fig.[T](a). 
The JC nonlincarity accurately describes the interaction 
of cavity QED systems both in the quantum (few photon) 
limit as well as in the semiclassical (many photon) limit 
(48j . In addition, it allows for an investigation of detuning 
effects between the cavity photons and the TLS present 
within each cavity. 

In our model, we consider either incoherent driving or 
coherent driving on the TLS, corresponding to two dis- 
tinct experimental situations. As a technical point, we 
introduce the Arnoldi time evolution method in order to 
solve numerically exact the time dependence of the non- 
equilibrium system. Specifically, wc evaluate the fluores- 
cence spectrum and the spectrum of the second-order cor- 
relation function for the TLS. This corresponds to an ex- 
perimental situation where the TLS has a non-negligible 
loss rate into the non-cavity or 'transverse' modes. We 
specifically consider emission into such non-cavity modes 
as they are both conceptually and experimentally easier 
to unambiguously distinguish from the emission stem- 
ming from cavity loss itself. 

The main goal of this work is to relate spectra mea- 
sured for the laser-driven and dissipative quantum sys- 
tem to the excitation energies of the non-dissipative 
quantum system. We show that the fluorescence spec- 
trum probes excitations between two adjacent particle 
number sectors, whereas the spectrum of the second- 
order correlation function probes density-density excita- 
tions within one particle number sector. The second- 
order correlation function of the transverse field emitted 
from cavity systems with a JC nonlincarity exhibits qual- 
itatively different properties when compared to systems 
with a Kerr nonlincarity, as the emitted photons are al- 
ways antibunched. 

The remainder of this paper is organized as follows. In 
Sec. U we introduce our model in detail, which describes 
a laser-driven and dissipative system of two coupled cav- 
ities with JC nonlinearity. In Sec. [II] we discuss the eval- 
uation of first- and second-order steady-state correlation 
functions and spectra using the quantum regression the- 
orem [54| . The Arnoldi time evolution, which is the nu- 
merical method used to calculate the steady-state corre- 
lation functions is introduced in Sec. IIIII Section IIVI is 
devoted to incoherent driving as a pump mechanism. In 
particular, single-particle excitation spectra and second- 
order correlation functions of the TLS situated in the sec- 
ond cavity are evaluated and discussed. Coherent driving 
is investigated in Sec.fV] Finally, we discuss and conclude 
our findings in Sec. IVII 



I. TWO LASER-DRIVEN AND DISSIPATIVE 
COUPLED CAVITIES 

In this paper, we investigate a laser-driven and dissi- 
pative system of two coupled cavities. A single cavity at 
site i is modeled by the JC Hamiltonian [4?} 

H- C = uj c a\ a ; + e af + g(a { a? + aj <7 t ~) , (1) 

where we set h = 1 as we do throughout this paper. The 
operator a\ (a i ) creates (destroys) a photon of frequency 
to c at lattice site i. The TLS is mathematically described 
by Pauli spin algebra, where af (a~ ) is the raising (low- 
ering) operator. The lower energy level will be denoted 
by \l) and the upper energy level by The energy spac- 
ing of the TLS is e, and g is the dipole coupling strength 
between photons and the TLS. In the subsequent analysis 
and calculations we use g as our unit of energy. When 
considering the JC model in the rotating frame of the 
photon frequency uj c , the only remaining parameter is 
the detuning A between the photons and the two-level 
excitations A = uj c — e. 

The total Hamiltonian of two coupled cavities is given 

by 
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H = -J(a\a 2 + a 1 al) + Y / H l JC . (2) 

i=l 

The coupling of the cavities arises due to the finite over- 
lap of their photonic wave functions, which determines 
the hopping strength J, see sketch in Fig.Q](a). 

In experiments cavity systems are driven by an ex- 
ternal laser source and are susceptible to photon loss. 
Therefore we have to investigate an open quantum sys- 
tem, whose d yna mics are recovered under certain approx- 
imations [5J-l56| by a master equation for the density 
operator p 



The superoperator C acts on the density operator p as 
tp = -i[H, p]+J2 I V /•,/'/•; - \{L\hiP)) • ( 4 ) 

The first term on the right hand side of Eq. (j4} accounts 
for the unitary time evolution and the second term for 
decoherence processes. The operator Lj is the Lindbla- 
dian corresponding to the jth decoherence channel with 
characteristic decoherence rate Tj . Square brackets indi- 
cate the commutator [a, b] = ab — b a and curly brackets 
the anticommutator {a, 6} = ab + ba. 

In our investigations, we consider two distinct kinds of 
laser driving corresponding to two different experimen- 
tal situations. The first is incoherent driving, which is 
sketched in Fig.[T](b) and takes advantage of a three level 
structure. An external laser drives the transition from \\.) 
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to \p) coherently. A part of the excitations in \p) decays to 
the metastable state |f) by spontaneous emission, which 
leads to occupation of the upper state of the TLS consid- 
ered in the JC Hamiltonian, see Eq. ([1]). This driving can 
be treated by an incoherent pump term, leading to the 
Lindbladian L\ = <r^ with effective pump rate Ti = T + , 
where we assume that the first cavity is driven by the 
external laser as indicated in Fig.[T](a). 

The second kind of driving we consider is coherent driv- 
ing, where the TLS is directly driven by a laser of fre- 
quency Ul, see Fig.[T](c). This leads to an additional 
Hamiltonian term y57f 



H 



drv 



Fe- iuJLt at + F*e luJLt c 



which must be added to H. When considering the driving 
strength F = T + to be real, which affects only a global 
phase, and writing the Hamiltonian in the rotating frame 
of the laser frequency ujl, H drv reduces to 



H 



drv 



(5) 



The impact of the rotating frame on the total Hamil- 
tonian H is that there exist two detunings, specifically, 
the detuning between the laser and the cavity photons 
8 u>c = uil — w c and the detuning between the laser and 
the energy spacing of the TLS S e = ujl — £• 

In cavity and circuit QED experiments, many loss 
channels have to be considered for photon decay. Among 
them are cavity losses and losses arising from decay of the 
TLSs into non-cavity modes. Our goal is not to describe 
all the details of a certain experiment, but to investi- 
gate the simplest possible model, that covers the essen- 
tial physics of laser-driven and dissipative cavity systems. 
Thus we consider only one loss channel coupling to the 
TLS of the second cavity, meaning we consider a Lind- 
bladian loss term L 2 = <J 2 with rate F2 = T~ . A sketch 
of the two coupled cavities with driving and dissipation 
is shown in Fig.[T](a). 



II. 



CORRELATION FUNCTIONS 



The purpose of this paper is to analyze spectra and 
correlation functions of photons emitted from the TLS 
of the second cavity. We investigate the TLS correlation 
functions, assuming that emission of the TLS into non- 
cavity modes can be differentiated from ordinary cavity 
losses. This process is best visualized as the emission of 
the TLS into modes that are transverse to the plane of 
the cavity confined photonic modes. 

All quantities arc evaluated in steady state p ss , which 
is obtained by solving Eq. ^ for dp/dt = 0, i.e., by 
solving Cp = 0. From the experimental viewpoint ob- 
servables evaluated in steady state, i. e., in the long time 
limit, are of direct relevance, as in most experiments 
there is enough time between turning on the lasers and 
performing the measurements for the system to reach 
equilibrium. 



In particular, we calculate the fluorescence spectrum, 
which is obtained from the Fourier transformation of the 
first-order correlation function in the steady state limit 



dt 

7^ 



"(*K"(o)>, 



(6) 



where (o- 2 (t)a 2 (0)) ss = \im T ^ 00 (<r 2 h (t + t) a 2 (r)). The 
correlation function can be conveniently evaluated using 
the quantum regression theorem For t > we have 



(a+(t) a 2 (0)} ss = Tr(<r+ e Lt a 2 p ss ) , 



(7) 



where exp[££] is the formal solution of the master equa- 
tion in Eq. (J3J. 

In addition to the fluorescence spectrum, we analyze 
the spectrum G(ui) of the second-order correlation func- 
tion g( 2 \t). The second-order correlation function for 
the TLS of the second cavity is defined as 



(<r+(0)<r+(t) a 2 (t)a 2 (0)) ss 
(a+(0)a 2 (0)) ss (a+(t)a 2 (t)), 



(8) 



We again employ the quantum regression theorem and 
obtain for t > 



9 m it) = 



_ Tt[<jJ a 2 e ct a 2 p ss a} 



(Tr[CT+ cr 2 p ss ]) 2 



(9) 



The second-order correlation function g( 2 > (t) provides in- 
formation about the statistics of the photons. Specifi- 
cally, it is the probability to detect two photons sepa- 
rated by a time delay t. In the case of a coherent light 
field g( 2 '(t) = 1 for all t indicating a completely random 
sequence of photon pulses leading to a Poissonian distri- 
bution of time intervals between the pulses. 

However, when g^ 2 '{0) < 1 the emitted light exhibits 
distinct quantum mechanical features, i.e., antibunching 
[HI, [HH . Antibunched light leads to a sub-Poissonian dis- 
tribution of time intervals between the emitted photons. 
The light emitted from the TLS considered here is per- 
fectly antibunched, due to the nature of a TLS, since 
(°2~) 2 = ( CT 2~) 2 — an d thus we have g^iO) = 0, see 
Eq. ([5J . We introduce the Fourier transformation of the 
second-order correlation function as 



dt 



e iut fl (2) (t) 



(10) 



which, as we show below, contains information about the 
density-density correlations. 



III. 



ARNOLDI TIME EVOLUTION 



In this section, we introduce the Arnoldi time evolu- 
tion for open quantum systems, to obtain the full quan- 
tum mechanical solution of the master equation given in 
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Eq. (|3|). To this end, we rewrite the master equation as 
a matrix equation of the form 



(11) 



where p = vec(p) is a column-ordered vector represen- 
tation of the density matrix and C is the superoperator 
in matrix form, which includes effects of both the uni- 
tary time evolution and the decoherence processes. Using 
tensor product properties [59l |60| the following identity 
holds 



vec(AXB) = (A(g>B T )X 



(12) 



Thus with judicious insertion of the identity operator 1 
into Eqs. © and (HJ) we have 



1-/51 = -i(Hpl-lpH) 
at 

+ Y. V ^ l j'P L ] - liL^pl + lpL^)] . (13) 



Applying the identity Eq. (|12p yields the wanted matrix 
representation of the master equation 



dt P 



-i[H® l-l®H T ]p 

]T ^[L. ® L* - ® 1 + 1 ^(Lt^mp 

(14) 



The superoperator matrix £ is of size N 2 x TV 2 and p is a 
vector of size A 2 , where N is the Hilbert space dimension 
of the two cavity system. 

The formal solution of Eq. (|TT|) is 



p(t) = e"p(0) 



(15) 



In practice, the observables and thus the density operator 
have to be evaluated at multiple times separated by At 
leading to the propagation prescription 



p(t + At) = e CAt p{t) 



(16) 



We are now left with the calculation of the exponen- 
tial exp[£Ai], which amounts to diagonalizing the large, 
sparse, complex valued and non-Hermitian matrix C. 
This can be achieved by the Arnoldi method [(3l| , which is 
an iterative technique to approximate the eigenvalues of 
such matrices. In the following, we briefly sketch how the 
Arnoldi method works. Details can be found in Refs. [62| 
and [6^ . This procedure is essentially based on construct- 
ing an orthonormal basis V of the Krylov subspacc 

lC m (£,p) = span{p, Cp, C 2 p, . . . C^p} . 

In the following we will refer to V as the Arnoldi ba- 
sis. The orthonormal basis V projects the superoperator 
matrix C onto an upper Hesscnberg form H 



V^CV = H 



(17) 



V 1 




V 



FIG. 2. (Color online) Illustration of the projection of the 
superoperator matrix £ onto an upper Hessenberg form H 
using the Arnoldi basis V . 



An illustration of this projection is shown in Fig. [5] The 
upper Hesscnberg form H is a small matrix of size m x m 
and can be diagonalized using standard methods 



HU = UD 



(18) 



where the matrix U contains the eigenvectors of H and 
the diagonal matrix D its ei genv alues, which are Ritz 
approximate eigenvalues of C [62||. Combining Eqs. (fl7|) 
and (|18|) we have 

C = VXJDXJ- X V^ . 

Since V is an orthonormal basis, i.e., V'V — 1, the 
formal solution given in Eq. (|16p can be expressed as 



p(t + At) = VU e DAt U- 1 V^p(t) . 

This expression is evaluated straightforwardly as the ma- 
trix D is diagonal. When choosing p(t) as initial vec- 
tor of the Krylov subspace, V'p(t) simplifies to e\ = 
(1, 0, . . .) T since all vectors of the Arnoldi basis V arc 
constructed to be orthogonal. With these considerations 
we have 



p(t + At) = VUe DAt U- 1 e 1 . 



(19) 



The Arnoldi method is an iterative procedure and thus 
a convergence criteria has to be imposed. The propa- 
gation Eq. (fl~9"|) can be regarded from a slightly different 
viewpoint as an expansion in the Arnoldi basis vectors 



p(t + At) = J2 



CiVi 



(20) 



where Hi is the ith column of V and Ci is the ith com- 
ponent of the vector c = U e DAt U~ x e\. A convenient 
convergence criteria for Krylov based time evolution 64 1 



\c m \ 2 <T, 



(21) 



since an extension of the Krylov basis has no implica- 
tions on the result within the order of the tolerance T. 
Importantly, the choice of both the time step At and the 
size of the Arnoldi basis m influences the convergence. In 
practice, the time step is chosen to be a trade off between 
the maximally reached time and the computer time (and 
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FIG. 3. (Color online) Fluorescence spectra S(uj) of the laser-driven and dissipative two-cavity system for hopping strength 
J/g = 0.2, loss rate F~ j g — 0.2 and detuning (a) A/g = and (b) A/g = 1, respectively. The driving strength T + is adjusted 
such that the steady-state particle density is n p = {0.25, 0.5, 0.75, 1.0}. In addition to the spectra, the light-gray (vertical) 
lines indicate single-particle excitations of the non-dissipative system from the N p = to N p — 1 particle sector of the Hilbert 
space, which we refer to as (0— s-l)-transitions. These excitations turn out to yield the main contributions to the spectra. The 
inset in panel (a) magnifies the two-peak structure of S(cj) at \u — uj c \ ~ —g and demonstrates the shift of the resonances with 
varying driving strength. 



memory) needed to perform one Arnoldi time step At, 
which of course depends on the Krylov basis size to. 



stochastically using quantum trajectories [5c 
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Alternatively to employing the Arnoldi time evolution, 
Eq. (jlip can be solved deterministically by conventional 
integration schemes based on Runge-Kutta like approxi- 
mations (for an application see e. g. Ref.l52l) or by Laplace 
transforming the differential equation and subsequently 
solving numerically a sparse system of equations for each 
frequency of interest [6a . |6(| . The latter approach is par- 
ticularly suitable for small systems, since it is a numer- 
ically exact approach and easy to apply. However, for 
medium sized or larger problems solving the sparse sys- 
tem of equations is not feasible any longer due to strin- 
gent memory requirements. Both conventional integra- 
tion schemes and the Arnoldi time evolution can be ap- 
plied to medium sized problems. Compared to conven- 
tional integration schemes, the Arnoldi time evolution 
method needs slightly more memory since, depending on 
the chosen timestep At, approximately 10 to 20 Krylov 
basis vectors have to be stored. The advantages of the 
Arnoldi time evolution, however, are that it is very stable 
in the long time limit and that a stringent convergence 
criteria [sec Eq. (|2"Tj) ] is available to control the accuracy 
of the evolved density matrix. Alternatively, to the deter- 
ministic approaches, the master equation can be solved 



In our calculations, we employ the Arnoldi time evolu- 
tion to evaluate the time dependent correlation functions 
from the full quantum mechanical master equation. The 
initial state for the time evolution is the steady-state den- 
sity matrix p ss appropriately multiplied by raising and 
lowering operators of the TLS, compare with Eqs. ([7]) and 
([9]). We evaluate the steady-state density matrix p ss by 
solving for the null space of the superoperator, i.e., by 
solving Cp = 0, which, in the case of our model, has a 
unique solution for any physically applicable set of pa- 
rameters. 



Technically, solving Cp = is achieved by a shift-and- 
invert Arnoldi method, where a small shift e regularizes C 
and the invert mode guarantees convergence toward the 
eigenvalue closest to e. For small enough shifts e, the 
shift-and-invert Arnoldi method converges toward the 
wanted eigenvalue of the superoperator £, and thus, 
the corresponding eigenvector is the steady-state density 
matrix in column order. 
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IV. 



INCOHERENT DRIVING 



(a) 



(b) 



In this section, we present results for the incoherently 
driven system of two coupled dissipative cavities. Specif- 
ically, we evaluate fluorescence spectra S(uj) and spectra 
G(w) of the second-order correlation functions g( 2 \i). 



A. Fluorescence spectra 

In order to understand quantum systems it is impor- 
tant to track properties gained from measurements on 
the realistic, dissipative system back to properties of the 
non-dissipative system. Here, we investigate fluorescence 
spectra of the dissipative quantum system and relate the 
resonance peaks of the spectra to excitations of the non- 
dissipative system. In Fig. [3] we show the fluorescence 
spectrum, i.e., the Fourier transform of the first-order 
correlation function, see Eq. ([6]) . 

We choose a small hopping strength J/g = 0.2, corre- 
sponding to the regime of strong correlations, since the 
infinitely large, non-dissipative coupled-cavity system ex- 
hibits a Mott to superffuid quantum phase transition at 
this parameter regime for small filling factors [U, 0, 0- 
[T2I . Ha , HH, EH ■ Thus most interesting physics is to be ex- 
pected in this parameter region. In addition, this param- 
eter regime is also of great interest for the two-cavity sys- 
tem, as can be seen from the particle-density diagram in 
the hopping strength and chemical potential plane eval- 
uated in Ref. 0. In this regime of the hopping strength, 
regions of considerable extent, for varying chemical po- 
tential, exist for arbitrary filling factors. The loss rate 
is set to T~/g = 0.2 and the strength of the driving 
laser T + is adjusted such that the steady-state density is 
rip = {0.25, 0.5, 0.75, 1.0}. As throughout this work we 
use the dipole coupling strength g as the unit of energy. 
In Fig. 0(a) the spectra are evaluated for zero detuning 
A/g = whereas in (b) a finite detuning of A/g = 1 is 
considered. 

We employ the Arnoldi time evolution described in 
Sec. IIIII to evaluate the time dependent correlation func- 
tions from the full quantum mechanical master equation. 
As a technical point, it is important to note that the time 
evolution can be performed only for positive time steps, 
since negative time steps would lead to numerical insta- 
bilities of the exponential in Eq. (|T9j) . However, negative 
times can be attained by observing that the first-order 
correlation function evaluated for a negative time —t is 
simply the complex conjugate of the one for positive time 
t, see Eq. ([7]). The Hilbert space of the two-cavity sys- 
tem is infinitely large, as each cavity can be occupied by 
infinitely many photons. For the numerical investigation 
we set the maximal number of photons per cavity to 8, 
which was found to result in convergence for the weak 
driving limit considered here. 

The analogy of the fluorescence spectra in condensed 
matter physics is the single-particle excitation spectrum, 
which probes the energy necessary to add (remove) a 



N p = 2 



N n = l 



N p = 
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FIG. 4. (Color online) Energy ladder for zero detuning A/g = 
of (a) a single cavity and of (b) two coupled cavities. The 
number of particles in the respective Hilbert space sector is 
N v . 



particle to (from) the system. The Hamiltonian of the 
non-dissipative two-cavity system can be arranged in 
disconnected sectors of a certain particle number N p . 
Due to the block diagonal form of the Hamiltonian the 
Schrodinger equation can be solved for each sector N p 
separately leading to 

H\^)=E^\^) , 

where E^ p are the cigenenergies and \ipu P ) are the cor- 
responding eigenvectors. The single-particle excitation 
spectrum probes excitations from the Hilbert space sec- 
tor with iVp-particles to (N p + l)-particles, correspond- 



In 



ing to the energy difference ui^S = E t 
the weak driving limit, we infer that the excitations from 
the zero-particle sector N p = to the one-particle sector 
N p = 1 are the most dominant. In the following, we will 
refer to these excitations as (0— >l)-transitions. The non- 
dissipative system can be solved analytically for these 
Hilbert space sectors, leading to 
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(22) 



In Fig. [4] we compare schematically the energy ladder 
for zero detuning A/g = of (a) a single JC cavity, 
obtained from the solution of Eq. (JT|) , and of (b) two cou- 
pled cavities. Most important for our analysis within 
the weak driving limit are the low energy excitations. For 
zero coupling J/g = the eigenenergies of the N p = 1 
particle number sector w c ± g are two-fold degenerate, 
see Eq. (j2"2"j) , and apparently correspond to the N p = I 
eigenenergies of a single JC cavity. A finite coupling 
strength J / g, however, removes the degeneracy and leads 
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FIG. 5. (Color online) Density plot of the fluorescence spectra S(to) for steady-state particle density n v = 0.5, top row, 
and n p = 1.0, bottom row. The loss rate is set to T~ / g = 0.2. The first two columns depict the density as function of the 
hopping strength J/g for two detuning parameters A/g — [column (a.*)] and A/g = 1 [column (b.*)]. The third column, 
labeled (c.*), shows the detuning-dependence of the fluorescence spectra for J/g — 0.2. In all panels dark dashed lines indicate 
(0— >l)-transitions of the non-dissipative system. 



to the four peak structure observed in Fig. [3] The ex- 
citations = — E® obtained from Eq. (|2"2"j) . corre- 
spond to the energy differences between these four levels 
and the ground state energy. In Fig.|3] these energy dif- 
ferences are indicated by light-gray (vertical) lines and 
match as expected the main four resonances, present in 
the spectrum of the dissipative quantum system. Impor- 
tantly, for increasing driving strength the resonances of 
the open quantum system are slightly shifted as demon- 
strated in the inset of Fig.[3](a). 

Fig. [5] shows density plots of the fluorescence spectra. 
The top and bottom rows depict results for steady-state 
particle densities n p = 0.5 and n p = 1, respectively. The 
first two columns show the spectra as function of the hop- 
ping strength J/g for two different detuning parameters 
A/g = [column (a.*)] and A/g = 1 [column (b.*)]. The 
third column, labeled (c.*), shows S(u) for J/g = 0.2 as 
a function of the detuning A. In all plots we set the loss 
rate to T~/g = 0.2. The (0— >l)-transitions of the non- 
dissipative system are indicated by dark dashed lines. In 
the case of weak driving the main contributions to the 
spectra of the dissipative two-cavity system correspond 
to these excitations. However, for increasing driving 
strength, comparing the top and bottom rows of Fig.[SJ 
excitations from higher particle number sectors become 
increasingly important. In the spectra shown here, the 
main contributions to the additional peaks stem from the 
(1— »2)-transitions. Interestingly, for increasing driving 
strength the weight of the (0— >l)-transitions stays almost 



constant, whereas the weight of the (1— >2)-transitions in- 
creases, in the regime considered here, since the N p = 1 
particle number sector has to be occupied before higher 
excitations can take place. When comparing Fig.[5](c.l) 
and (c.2) it can be observed that the detuning mainly 
affects the excitation energies of the (0— >T) -transitions. 
The excitation energies of the (1— >■ 2) -transitions remain 
almost unaltered, only the intensities of S(u) are slightly 
modified. 

The general conclusion of Fig. [5] is that all poles u^S 
resulting from excitations from the sectors with N p to 
N p + 1 particles contribute to the fluorescence spectra 
and the strength of the poles depend on the occupation 
probability of states associated with the transition. 



B. Approximate semianalytic spectrum 

To gain more insight into the excitations we approxi- 
mately evaluate the spectrum in a semianalytic way from 
the steady-state density matrix p ss . To this end, we as- 
sume that each resonance u a = i 



N„ 



where a is a collective index containing fi, v and N p , of 
the spectrum can be described by a Lorentzian, which is 
of the form 



S a (u}) = l c 



(r/2) 2 



{u-u a y + {Y/2f 
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where l a is the strength of the individual peak and T is 
its half-width. We evaluate the half-width V by 



exact 

semianalytic 



r = r - 



ph 



where T~ is the loss rate and n\ (nf ) is the TLS (pho- 
ton) steady-state density of the second cavity, where the 
emission spectrum is measured. Next, we have to find 
an approximate prescription to evaluate the individual 
strength l a of each pole u) a . Motivated by the Green's 
function theory for closed quantum systems [69j . we as- 
sume that the contribution of each pole is proportional 
to 



JVp + l| 



WI\W')(W'\<tZp..M' +1 )\ 



compare with the first-order correlation function given in 
Eq. (J7J . Following these considerations a Fourier trans- 
form of the approximated correlation function would 
yield a resonance peak located at the assumed energy 
u> a , which once again proves our choice of l a to be rea- 
sonable. The full spectrum is constructed by summing 
over all individual Lorentzians 



Results obtained from this semianalytic treatment are 
shown along with the exact numerical results in Fig. [6] for 
(a) zero detuning A/g = and (b) detuning A/g = 1. 
The hopping strength is J/g = 0.2, the loss rate is 
r~ / g = 0.2 and the driving strength T + is adjusted such 
that the steady-state particle density is n p = 0.5. We 
normalize the semianalytic spectrum such that the area 
under the curve coincides with the area under the ex- 
act numerical spectrum. Overall the spectrum evaluated 
with this semianalytic technique fits the exact numer- 
ically evaluated spectrum rather well. In addition, to 
the four major resonances, peaks corresponding to exci- 
tations from Hilbert space sectors with N p > can be 
observed. 



C. Second-order correlation functions 

Returning to our full numerical solution of the master 
equation, we are able to evaluate the time-dependence 
of the second-order correlation function g^ 2 \t) and its 
spectrum G(w). We evaluate g^it) for the TLS of the 
second cavity, since the transverse field emission can be 
conveniently separated from the emitted cavity photon 
field in experiments. The second-order correlation func- 
tion g^(t), see Eq. ([5]), is the probability of detecting two 
photons separated by a time delay t. As required for the 
second-order correlation function of a TLS g^(0) = 0, 
corresponding to a completely quantum mechanical, an- 
tibunched light source. This is in stark contrast to re- 
sults obtained from a Hubbard model with Kerr nonlin- 
carity, where a transition from bunched to antibunched 




FIG. 6. (Color online) Comparison between the exact spec- 
trum, red, solid line, and the semianalytic spectrum, blue, 
dashed line. The hopping strength is J/g = 0.2, the loss rate 
is r~ / g = 0.2 and the driving strength L + is adjusted such 
that the steady-state particle density is n p — 0.5. Spectra 
are evaluated for detuning (a) A/g = and (b) A/g = 1. 
Light-gray (vertical) lines indicate (0— >l)-transitions of the 
non-dissipative system. 



light statistics can be observed as a function of the non- 
linearity strength f38j - l45| . The second-order correlation 
function g^ 2 \t) is an oscillating function, which initially 
increases and eventually approaches its asymptotic value 
one, meaning that the two emitted photons are uncorre- 
cted for large time delay t. Such asymptotic behavior re- 
sults in a pronounced peak in G(oj) at w = 0. In Fig.[7]wc 
show the spectrum G(oS) for zero detuning A/g = 0, loss 
rate T~/g = 0.2, and driving strength r + adjusted such 
that the steady-state particle density is n p = {0.5, 1.0}. 
Results are evaluated for hopping strength (a) J/g — 0.2, 
(b) J/g = 0.4, and (c) J/g = 0.6. 

Due to the similarity of the second-order correlation 
function, given in Eq. ([8"j). with the density-density cor- 
relation function (^2(^)^2 where = cr^cr^, one is 
prompted to analyze the results of the second-order corre- 
lation function in terms of the density-density correlation 
in the non-dissipative system. The density-density cor- 
relation correspond to excitations within a Hilbert space 
sector of fixed particle number N p and therefore probes 
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FIG. 7. (Color online) Spectrum G(cj) of the second-order 
correlation function, for zero detuning A/g — 0, loss rate 
r~ jg = 0.2 and hopping strength (a) J/g = 0.2, (b) J/g — 
0.4, and (c) J/g = 0.6, respectively. The driving strength 
r + is adjusted such that the steady-state particle density is 
n p = {0.5, 1.0}. Light-gray (vertical) lines indicate (1—^1)- 
transitions within the N p = 1 sector of the non-dissipative 
two cavity system. 



(N p — >-/Vp)-transitions. The corresponding resonances are 

-N p „N„ 



located at oj^S = E p p — E v p . Light-gray (vertical) lines 
in Fig. [7] indicate the position of the poles originating 
from (1— >-l)-transitions of the non-dissipative two-cavity 
system Q x = Ej t — El as obtained from Eq. (|2"2"j) . The 
spectrum G(ui) exhibits resonances at the poles u)^ along 



with a background that stems from density-density ex- 
citations ui^S for the Hilbert space sectors with larger 
particle number N p > 1. 



V. COHERENT DRIVING 

Here, we investigate coherent driving on the TLS mod- 
eled with the additional Hamiltonian term introduced in 
Eq. ([5} . In this section, we restrict our analysis to zero 
detuning between the TLS and the cavity photons, i.e., 
uj c = e and thus A/g = 0. From this restriction it fol- 
lows that the detuning <5 e between the driving laser and 
the TLS is equal to the detuning 8 Uc between the driving 
laser and the cavity photons. In the following, we denote 
this coherent driving detuning as 8 C = 5 f = S Ua . 

The coherently driven system is populated only when 
the driving laser is in resonance with a transition of the 
two-cavity system. For vanishing loss rate and very weak 
driving only the N p = 1 sector of the initially empty 
two-cavity system can be occupied, since all transitions 
involving higher particle numbers are detuned from the 
(0— >l)-transitions and are thus prohibited. However, the 
dissipation is responsible for an energy level broaden- 
ing as compared to the non-dissipative system shown in 
Fig-Hl(b) and therefore has a drastic impact on the occu- 
pation. Following these considerations, the energy spec- 
trum of the coherently driven two-cavity system can be 
probed by measuring the occupation of the TLS in the 
second cavity 77*2 cLS cl function of the coherent driving 
laser detuning 5 C , which is accessible in experiments by 
measuring the output field jU H(| . 

Numerical results for the steady-state occupation nf 
as a function of the coherent driving detuning S c are de- 
picted in Fig. [51 main panel (a), for hopping strength 
J/g = 0.2, and rates {T+ / g = 0.04, Y~/g = 0.5}, 
{T+/g = 0.04, T-/g = 0.1}, and {T+/g = 0.01, V'/g = 
0.1}, see Eq. ©• A magnification of the positive detuning 
response for the respective rates are shown in panels (b)- 
(d). For large loss rate Y~ = 0.5 and driving r + = 0.04, 
see panel (b), the two peak structure expected for (0—^1)- 
transitions cannot be resolved due to the level broaden- 
ing and only one peak at \uj — oj c \ w g is observed in the 
occupation spectrum. When decreasing the loss rate to 
T~ = 0.1, and keeping the driving fixed at T + = 0.04, 
see panel (c), the expected four peak structure emerges. 
However, apart from the (0— >l)-transitions, which are in- 
dicated by solid light-gray (vertical) lines, additional res- 
onances at \uj — u) c \ w 0.7g and \lo — uj c \ w g are present. 
These resonances are identical to half of the eigenenergies 
of the N p = 2 particle number sector and originate 
from subsequent (0— >T)- and (1— »2)-transitions. The en- 
ergies are indicated by dashed light-gray (vertical) 
lines in panel (c) . In the non-dissipative case these tran- 
sitions are not allowed because each (0— »l)-transition is 
detuned from the (1— >-2)-transitions. However, for large 
enough driving strength and for finite loss rate, i.e., fi- 
nite level broadening, these transitions become possible. 
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FIG. 8. (Color online) Steady-state occupation n% of the TLS in the second cavity, main panel (a), as a function of the coherent 
driving laser detuning 8 C for hopping strength J/g = 0.2 and distinct values of the driving strength F + j g and loss rate Y~ /g, see 
legend. Panels (b)-(d) magnify the occupation n% for the positive detuning response of the individual rates. Solid, light-gray 
(vertical) lines indicate the eigenenergies El of the N p — 1 particle number sector of the non-dissipative two cavity system, 



whereas dashed, light-gray (vertical) lines in panel (c) 
number sector; see main text for details. 



correspond to half of the eigenenergies of the N p = 2 particle 



The additional peaks, however, disappear when reducing 
the driving strength to T + = 0.01, as shown in panel (d), 
where the only remaining peaks stem from the (0— >1)- 
transitions. 

The fluorescence spectrum S(uj) of a coherently driven 
system is non-zero only when the driving laser is in reso- 
nance with the two-cavity system, meaning that the driv- 
ing laser is able to populate the two-cavity system. Im- 
portantly, a resonance peak in the fluorescence spectrum 
is observed at the specific frequency imprinted from the 
detuning of the coherent driving laser S c . The only ex- 
ception takes place for additional transitions that emerge 
due to the level broadening. 



VI. CONCLUSIONS 

In the present paper, we investigated a laser-driven and 
dissipative system of two coupled cavities with Jaynes- 
Cummings nonlincarity. In particular, we presented 
and discussed the fluorescence spectrum S(uj) and the 
spectrum G(oS) of the second-order correlation function 
(t) evaluated in steady state, which is of utmost rel- 
evance for experiments. All quantities are obtained by 
solving the time dependence of the full quantum me- 
chanical master equation numerically by means of the 
Arnoldi time evolution method introduced here to treat 
open quantum systems. 

Each cavity of the investigated system is described by 
the Jaynes-Cummings model and thus confines photons 



and contains a two-level system. The two-level system 
couples to the photons, which gives rise to the nonlinear- 
ity The coupling between the two cavities results from 
a finite overlap of respective photonic wave functions. 
In order to incorporate dissipation and driving we em- 
ployed a master equation in Lindblad form. Specifically, 
we consider two experimentally distinct kinds of driving 
on the two-level system of the first cavity: (i) incoherent 
driving, which takes advantage of a three level structure, 
that is modeled by a Lindbladian pump term in the mas- 
ter equation approach and (ii) coherent driving which is 
modeled by an additional Hamiltonian term. The dissi- 
pation is taken into account by a Lindbladian loss term 
coupling to the two-level system of the second cavity. 

In this work, we focused on relating the resonance 
peaks of the measured spectra to the energy spectra of 
the non-dissipative two-cavity system. In the case of 
incoherent driving, resonance peaks in the fluorescence 
spectra could be related to excitations from the N p - to 
(N p + l)-particle sector of the non-dissipative two-cavity 
system. In the weak driving limit the main contributions 
emerge from excitations corresponding to the N p = 
to N p = 1 transitions. For increasing driving strength, 
however, excitations from higher particle number sectors, 
such as Np = 1 to N p = 2 transitions become more pro- 
nounced. In order to gain insight into the measured spec- 
tra, we also presented an approximate prescription to 
calculate the spectrum semianalytically from the steady- 
state density matrix and from the excitation energies of 
the non-dissipative system. Finally, we demonstrated 
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that the spectrum G(u>) of the second-order correlation 
function probes excitations within a certain particle num- 
ber sector. Therefore it is related to the density-density 
correlations. For coherent driving the driving laser has 
to be tuned to be in resonance with a transition of the 
two-cavity system. This is necessary to populate the sys- 
tem. Thus we mapped out the energy spectrum of the 
non-dissipative two-cavity system by measuring the occu- 
pation of the two-level system in the second cavity. In the 
case of resonance, the spectrum yields a response at the 
detuning between the driving laser and the system, oth- 
erwise the spectrum is virtually zero. Additional peaks 
may occur due to the finite level broadening arising from 
the photon loss term. In summary, for coherent driv- 
ing, the occupation number provides a simple approach 



to obtain the excitation energies of the non-dissipative 
two-cavity system for small filling factors, whereas for 
incoherent driving, the spectra yield the full information 
about the excitation energies for arbitrary filling factors. 
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